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Abstract 

We study a class of random matrices that appear in several communication and signal processing 
applications, and whose asymptotic eigenvalue distribution is closely related to the reconstruction error of 
an irregularly sampled bandlimited signal. We focus on the case where the random variables characterizing 
these matrices are d-dimensional vectors, independent, and quasi-equally spaced, i.e., they have an 
arbitrary distribution and their averages are vertices of a cJ-dimensional grid. Although a closed form 
expression of the eigenvalue distribution is still unknown, under these conditions we are able (f) to derive 
the distribution moments as the matrix size grows to infinity, while its aspect ratio is kept constant, and 
(if) to show that the eigenvalue distribution tends to the Marcenko-Pastur law as d — > oo. These results 
can find application in several fields, as an example we show how they can be used for the estimation 
of the mean square error provided by linear reconstruction techniques. 

EDICS: DSP-RECO Signal reconstruction, DSP-SAMP Sampling, SPC-PERF Performance analysis 
and bounds. 
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I. Introduction 

Consider the class of random matrices of size (2M + 1) x r, with entries given by 



1 



V2M + 1 



-]2iTMx r 



(1) 



The generic element of G can be written as: G^ 



V2M+1 

where x g are independent random variables characterized by a probability density function (pdf) f x (z), 
with < z < 1. These matrices are Vandermonde matrices with complex exponential entries; they appear 
in many signal/image processing applications and have been studied in a number of recent works, (see 
e.g., [l]-[8]). More specifically, in the field of signal processing for sensor networks, [1], [2] studied 
the performance of linear reconstruction techniques for physical fields irregularly sampled by sensors. In 
such scenario, the random variables x q in (Q]) represent the coordinates of the sensor nodes. The work 
in [3] addressed the case where these coordinates are uniformly distributed and subject to an unknown 
jitter. In the field of communications, the study in [8] presented a number of applications where these 
matrices appear, which range from multiuser MIMO systems to multifold scattering. 

In spite of their numerous applications, few results are known for the Vandermonde matrices in (Q]). In 
particular, a closed form expression for the eigenvalue distribution of the Hermitian Toeplitz matrix GG^, 
as well as its asymptotic behavior, would be of great interest. As an example, in [1], [2], [6], it has been 
observed that the performance of linear techniques for reconstructing a signal from a set of irregularly- 
spaced samples with known coordinates is a function of the asymptotic eigenvalue distribution of GG^. 
The asymptotic eigenvalue distribution of GG^ is defined as the distribution of its eigenvalues, in the 
limit of M and r growing to infinity while their ratio is kept constant. Unfortunately, such distribution 
is still unknown. 

In this work, we consider a general formulation which extends the model in ([T]) to the <i-dimensional 
domain. We study the properties of random matrices of size (2M + l) d x r and entries given by 



1 



V(2M + lF 



-j27r£ T x, 



(2) 



where the vectors x g = [x q ±, . . . , x q d] T have independent entries, characterized by the pdf f x m (z), 
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<7 = 0,... , r — 1, m = 1, . . . , d, and d is the number of dimensions. The invertible function 

d 

v(£) = J2(2M + l) m - 1 e m (3) 

m=l 

maps the vector of integers £ = [£x, . . . , £d] T , £ m = —M, ... ,M onto a scalar index, i.e., the row index 
of the matrix G^. Notice that, when d = 1, Gd reduces to (fl}. 

For the matrix model in Q, we study the interesting case where are independent, quasi-equally 
spaced random variables in the (i-dimensional hypercube [0, l) d . In other words, we assume that the 
averages of x 9 are the vertices of a (i-dimensional grid in [0, l) d . This is often the case arising in 
measurement systems affected by jitter, or in sensor network deployments where the sensors sampling 
the physical field can only be roughly placed at equally spaced positions, due to terrain conditions and 
deployment practicality [9]. Note that the distribution of the random variables x can be of any kind, the 
only assumption we make is on their averages being equally spaced. Since an analytic expression of the 
eigenvalue distribution of G^G^ is unknown, we derive a closed form expression for its moments. This 
enables us to show that, as d — > oo, the eigenvalue distribution tends to the Marcenko-Pastur law [12]. 
At the end of the paper, we present some numerical results and applications where the moments and the 
asymptotic approximation to the eigenvalue distribution of G^G^ can be of great use. 



As a first step, we briefly review previous results on the G^ matrices. In a one-dimensional domain 
(d = 1), the work in [1] considered an irregularly sampled bandlimited signal, which is reconstructed 
using linear techniques and assuming the samples coordinates to be known. The performance of the 
reconstruction system was derived as a function of the eigenvalue distribution f\(l,(3,z) of the matrix 
Ti = (3GiG\, where (3 is the aspect raticQ of Gi [1], [2]. An explicit expression of the moments 



was attained in [4], [5], for the specific case where x q are uniformly distributed in [0, 1). Also, in the 
case where x q are independent, quasi-equally spaced random variables, the analytic expression of the 
second moment of the eigenvalue distribution of T, i.e., E[Af J, was obtained in [3]. Then, in [7] the 
moments f\(l,f3,z) were derived for an arbitrary distribution f x (z). 

In [4], [5], the (i-dimensional model (O was also investigated. There, the properties of the random 
matrices G^ were studied in the case where the vectors x g = [x q \, . . . ,x q( i] have independent entries, 

'The aspect ratio of G is the ratio between the number of rows and the number of columns of the matrix 



II. Previous Results and Problem Formulation 
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uniformly distributed in the hypercube [0, l) d . Under such assumptions, and for given d and aspect ratio 
[5, an analytic expression of the moments of f\(d,(3, z) was derived and it was shown that, as d — > oo, 
f\(d,(3,z) tends to the Marcenko-Pastur law [12], i.e., 

r f <a r \ t fa \ V( c i -z){z-c 2 ) 
hm fx{d,(3,z) = f M p{P,z) = — - 

d-^oo ZTTZp 

where ci,c 2 = (1 ± \//3) 2 , < (3 < 1, c 2 < x < c\. 

The following sections detail the problem addressed in this work and introduce some useful notations. 

A. The quasi-equally spaced multidimensional model 

We consider the matrix class in (J2J) and assume that the vectors x are independent, quasi-equally 
spaced random variables in the d-dimensional hypercube [0, l) d , i.e., the averages of x are the vertices 
of a ci-dimensional grid in [0, l) d . 

We define p as the number of vertices per dimension, thus the total number of vertices is r = p d . We 
denote the coordinate of a generic vertex of the grid by the vector q/p £ [0, l) d , where q = [qi, . . . , qd] T , 
is an integer vector and q m = 0, . . . , p — 1. For notation simplicity and in analogy with (O, we identify 
the vertex with coordinate q/ p by the scalar index 

d 

Mq) = P" 1 ' 1 ^ < 4 > 

m=l 

Notice that < /i(q) < r — 1 is an invertible function and allows us to write 

q x M(q) 

X M(q) = - H 

p p 

where the average 

E[X ^) ] = ~p + Yp 

is the coordinate of the sample identified by the scalar label p(q) and 1 is the all ones vector. Furthermore, 
we assume that the entries of the vectors x^( q ) are i.i.d. with pdf fx{z) which does not depend on r, M, 
or q. By using this notation, the entries of are then given by 

1 p -j2^£ T x 

v/(2M + l) d 

while the aspect ratio is 



(GdW),Mq) = 7= . ^ e" j2 ^ ^ (5) 



^_ (2M + l) d _ / 2M + l \ d 
r VP/ 



The Hermitian Toeplitz matrix = /3GdG d is defined as 

q 
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where ^ represents a ci-dimensional sum over all vectors q such that q m = 0, . . . , p — 1, m = 1, . . . , d. 

Our goals are (i) to derive the analytic expression of the moments of f\(d,(3,z) with quasi-equally 
spaced vectors (Section |III]), and (ii) to show that as d — > oo, f\(d,/3,z) tends to the Marcenko- 

Pastur law (Section HVl) . 

III. Closed form expression of the moments of the asymptotic eigenvalue pdf 

Following the approach adopted in [13], [14], in the limit for M and r growing to infinity with constant 
aspect ratio (3 and dimension d, we compute the closed form expression of K[X^ J, which can be obtained 
from the powers of as [15], 

Tr{ E [T p d ] } 

E[A^ J = lim — — rr (8) 

L d >P s M,r->+oc (2M + l) d 

In <(8]) the symbol Tr identifies the matrix trace operator, and the average EH is computed over the set 

x 

of random variables X = {x , . . . , x r _i}. An efficient method to compute © exploits set partitioning. 
Indeed, note that the power T p d is the matrix product of p copies of T^. This operation yields exponential 
terms, whose exponents are given by a sum of p terms of the form x ^(q.)(^i — %+i]) ( see a ^ so in 
Appendix |A]). The average of this sum depends on the number of distinct vectors q^, and all possible 
cases can be described as partitions of the set V = {1,... ,p}- In particular, the case where in the 
set {qi, . . . ,q p } there are 1 < k < p distinct vectors, corresponds to a partition of V in k subsets. It 
follows that a fundamental step to calculate ([8]) is the computation of all possible partitions of set V. 
Before proceeding further in our analysis, we therefore introduce some useful definitions related to set 
partitioning. 

A. Definitions 

Let the integer p denote the moment order and let the vector /x = [p±, . . . , p p ] be a possible combination 
of p integers. In our specific case, each entry of the vector /x is given by the expression in (01), i.e., 
Pi = p(q_i) and, thus, can range between and r — 1. 

We define: 

• the scalar integer 1 < k(p) < p as the number of distinct entries of the vector /x; 

• 'y(fJ-) as the vector of integers, of length k(fi), whose entries 7j (/x), j = 1, . . . , k(fi), are the entries 
of fx without repetitions, in order of appearance within /x; 

• Vj(fM) as the set of indices of the entries of /x with value 7j(/x), j = 1, . . . , /c(/x); 
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• the vector wQu) = [ooi(fx), . . . , to p (fi)] such that, for any given j = 1, . . . , k(fx), we have = j 
if i G Vj{n), i = l,..., p. 

Example 1: Let fx = [1,5,2,8,5,3,2], then k{fx) = 5 since the entries of fx take 5 distinct 
values (i.e., {1,5,2,8,3}). Such values, taken in order of appearance in fx form the vector 
7(/x) = [1,5,2,8,3]. The value 71 = 1 appears at position 1 in fx, therefore Pi(fx) = {1}. 
The value 72 = 5 appears at positions 2 and 5 in ft, therefore T^U*) = {2,5}. Similarly 
Ps(fi) = {3,7}, Pi(fx) = {4}, and P^fx) = {6}. By using the sets Pj we build the vector, 
cj(/i). For each j = 1, . . . , k we assign the value j to every u>i such that i G Vj. For example, 
^2 = CJ5 = 2 since the integers 2 and 5 are in Vi- In conclusion cj(/i) = [1, 2, 3, 4, 2, 5, 3]. 

Furthermore, we define: 

• ftp as the set of partitions of V; 

p 

• VLp k as the set of partitions of V in k subsets, 1 < k < p, with U f2 Pi fc = f2 p . 

fc=i 

Note that: (i) the cardinality of J7 P , denoted by = |n p |, is the p-th Bell number [16] and (ii) the 

cardinality of £l Pt k, denoted by S(p,k) = |fi p ,fc|, is a Stirling number of the second kind [17]. 
From the above definitions, it follows that: 

1) the vector fx induces a partition of the set V which is identified by the subsets Vj{n). These subsets 
have the following properties 

Mm) 

Even though the partition identified by fx is often represented as {Pi, . . . , Pk(n)}, by its definition, 
an equivalent representation of such partition is given by the vector cj(/x). Therefore, from now on 
we will refer to cj(/x) as a partition of the p element set P induced by fx (for simplicity, however, 
often we will not explicit the dependency of u on fx); 

2) fc(aj) = k(fx), since the entries of u take all possible values in the set {1, . . . , k(fx)}; 

3) P j {u)=P j {fx), for j = l,..., k{fi). 

At last, we define M.(u>) as the set of fx inducing the same partition cj of P. 
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Example 2: Let r = 3 and p = 3. Since /x = [/zi , . . . , fj, p ] and fa = 0, . . . , r — 1, i = 1, . . . ,p, 
we have r p = 27 possible vectors n, namely, {[0, 0, 0], [0, 0, 1], . . . , [2, 2, 1], [2, 2, 2]}. Each 
fi identifies a partition cj G ^3.fc, with /c = 1, ... ,3, as described in Example 1. The sets 
of partitions 3)fc , are given by f2 3j i = {[1,1,1]}, Q 3:2 = {[1, 1, 2], [1, 2, 1], [1, 2, 2]}, and 
^3,3 = {[1; 2, 3]}, and have cardinality 5(3, 1) = 1, 5(3, 2) = 3 and 5(3, 3) = 1, respectively. 
The set of vectors /i identifying the partition u> = [1,1,1], i.e., 7W([1,1,1]), is given by: 
M([l, 1,1]) ={[0,0,0], [1,1,1], [2, 2, 2]}. Similarly, 

M([l, 1, 2]) = {[0, 0, 1], [0, 0, 2], [1, 1, 0], [1, 1, 2], [2, 2, 0], [2, 2, 1]} 

M([l, 2, 1]) = {[0, 1, 0], [0, 2, 0], [1, 0, 1], [1, 2, 1], [2, 0, 2], [2, 1, 2]} 

M([l, 2, 2]) = {[0, 1, 1], [0, 2, 2], [1, 0, 0], [1, 2, 2], [2, 0, 0], [2, 1, 1]} 

M([l, 2, 3]) = {[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]} 



B. Closed form expression of E[A^ g] 

By using the definitions in Section IIII-AI and by applying set partitioning to ([8]), we can state the first 
main result of this work: 

Theorem 3.1: Let be a (2M + l) d x (2M + l) d Hermitian random matrix as defined in (0, where 
the properties of the random vectors x^( q ) are described in Section III- A I Then, for any given f3 and d, 
the p-th moment of the asymptotic eigenvalue distribution of is given by: 



p k 



where 



k=l h=l ueHp.k Lj'eQ k:h 



u(u') = (-l) k - h \{{\V r (u')\-l)\ 

j'=l 



(9) 



(10) 



v\u, u 



h = 1 



/ f[C ^n^Wjiu)) dy 

k h / _ \ 

/ n C (-W l/ S( w ))[[ & £ «»i'H dy Kki (11) 



J|<5 D (twj-(a;)) dy 



h = k 



p j=i 
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TABLE I 

Partition sets Q. n , m for n = 1, 2, 3, and 1 < m < n. Each partition is represented through its associated 

VECTOR OJ AND THE VALUE OF u(w) 



ui, u(u>) 


m = 1 


m = 2 


m = 3 


n = 1 


[1], 1 






n = 2 


[1,1], -1 


[1,2], 1 




n = 3 


[1,1,1], 2 


[1,1,2], -1 
[1,2,1], -1 
[1,2,2], -1 


[1,2,3], 1 



and f(oj,cj') = 1 for k = 1. In (fTTI) . we denned W p as the p-dimensional hypercube [— 1/2, l/2) p , 
C(s) = E[e sz ] as the characteristic function of x, Sd(-) as the Dirac's delta, and 

X 

Wj{u)= Vi-V[i+1] 

yi e R, i = 1, . . . ,p, and j = 1, . . . , 

Proof: The proof can be found in Appendix |A] ■ 

With the aim to give an intuitive explanation of the above expressions, note that the right hand side 
of (© counts all possible partitions of the set V = {1, . . . ,p}, C(s) in (fTTI ) accounts for the generic 
distribution of the variables x, and the quantity Wj(u) represents the indices pairing that appears in the 
exponent of the generic entry of the power T^. 

To further clarify the moments computation, Table |TJ reports an example of partition sets Q njm for 
n = 1, ... ,3 and 1 < m < n, while Example 3 shows the computation of the second moment of the 
eigenvalue distribution. 
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Example 3: We compute the analytic expression of E[A^ J. Using (O, we get 



* = EE^ E E <"'>(»,"') d 

k=l h=l Lj£fl 2tk u}'eQ k:h 

By expanding this expression and using Table H we obtain 

E[A^] = f3v([l, 1], [l]) d - /3v([l, 2], [1, l]) d + v([l, 2], [1, 2]) d 

We notice that, for k = 1, u([l, 1], [1]) = 1. The term v([l, 2], [1, 2]) refers instead to the case 
k = h = 2, and it is given by 



r 2 

w([l,2],[l,2]) = / l[S D ( Wj ([l,2]))dy 



with ioi([l, 2]) = yi — U2 and u^Ql, 2]) = 2/2 — 2/1- It follows that 

v([l, 2], [1, 2]) = / - y 2 )M2/2 - dy = 1 

J'Hi 



Finally, 



([1,2], [1,1]) = I nc(-j2vr/3 1 /^.([i, 2 ])) dy 



Thus, we write 



E[A^] = l + /3-/3 



n 2 



n-2 



C[-]2TTf3 l l d {y x -y 2 



dy 



C (-j2vr/? 1 / d (y 1 - y 2 )) 



dy 



IV. Convergence to the Marcenko-Pastur distribution 

In this section we show that the asymptotic eigenvalue distribution of the matrix tends to the 
Marcenko-Pastur law [12], as d —> 00. This is equivalent to prove that, as d — > 00, the p-th moment of 
\d,p tends to the p-th moment of the Marcenko-Pastur distribution with parameter j3, for every p > 1. 

Theorem 4.1: Let be a (2M + l) d x (2M + l) d Hermitian random matrix as defined in ([7]), where 
the properties of the random vectors x^( q ) are described in Section ITl-AI Let E[A^ a] be the p-th moment 
of the asymptotic eigenvalue distribution of T^, given by Theorem 13.11 Then, for any given (3, 

M^]=E[^j] = j2P p - k N{p,k) (12) 
00 k=i 



June 23, 2008 



DRAFT 



10 



where N(j>, k) are the Narayana numbers [18], [19] and ElA^ p] are the Narayana polynomials, i.e., the 
moments of the Marcenko-Pastur distribution [12]. 
Proof: 

We first look at the expression of the p-th asymptotic moment and observe that, for h = k, the 
contribution of the term in the right hand side of Q reduces to 

E E «(«>(« V)* (13) 

k=l LJ£fl P ,k u'£Qk,k 

The cardinality of Q^k is ^(^j k) = 1 and 0^. = {[1, . . . , k]}. Thus, we only consider u>' = [1, . . . , fe]. 
Moreover, using ( fTOl ) we have u([l, ...,&]) = 1 since each subset Vj>([l, ■ ■ ■ , k]) has cardinality 1, 
f = 1, . . . , k. Therefore, the term in (fT3l) becomes 

fe=l uiSOp fc 

Using (fTTb with /i = k, we have: 

v(u,[l,...,k})= f JlM^C"))^ = < 14 > 

■'W. i=l 

Hence, the contribution to the p-th moment reduces to 

f>P- fc J] (15) 

fc=l OJgilp.fc 

In [4], [5] it is shown that, as d — ► oo, (031 ) tends to the Narayana polynomial of order p. It follows that, 
in order to prove the theorem, it is enough to show that for h < k the contribution of the term in the 
right hand side of ©, to the expression of the p-th asymptotic moment, vanishes as d — > oo. In practice 
we have to show that, for each u 6 Q Pj k and u>' G £lk,h> with h < k, 

lim v{u,u') d = 

d— >oo 

or, equivalently, that |v(u>,u/)| < 1. 
We first notice that for 1 < h < k 
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Moreover, we have: 

C'(-j2^ 1/ S(cj) 



exp(-j2vr/3 1 %.( Ci ,)^/-( z )dz 

-oo ^ ' 

/ + 0O , . 

ew(-&*P 1/d Wj(u)*) fx{*) 
-oo 

/+oo 
h(z)dz = l 
-oo 



(a) f+oo 
< 



(17) 



The equality (a) arises if the condition Wj(u) = is always verified, otherwise, if Wj{u>) ^ 0, 
\C (-i2irl3 l l d Wj {u)) \ < 1. 

Next, we make the following observations: (/) since we consider partitions u' of the form {1, . . . , k} 
in h subsets with h < k, then at least one of the sets Vj'(uj') has cardinality \Vj>(uj')\ > 1; (ii) the term 

H ( 

gives a non-zero contribution to the integral in (fl6l ) only when Yli'EV-,(u>') w i'( u; ) = 0- Hence, if 
|"Pj/(u/)| > 1 for some j', then some Wi>(u') ^ will provide a non-zero contribution to the integral 
in (fl6l ). In this case, we can write 

. k h I \ 

kwvji < / nl^H^s-H)! e m«)U 

Jn *j=i i'=i Vi'ePj/Cw') / 

< / IIm E ^'( cj )( dy^ 1 

which proves the claim. 

When ft = 1, again, there is a measurable subset of TL P for which Wj(u) 7^ 0, hence, 



(18) 



Kw,a/)| < ! n|c(-j2vr/3 1 /S-M 
i.e., the strict inequality holds. 



dy < 1 



In Figure [Q we show the empirical eigenvalue distribution of the matrix for (5 = 0.55, d = 1,2, 3, 
and x uniformly distributed in [0, 1]. The empirical distribution is compared to the Marcenko-Pastur 
distribution (solid line). We observe that as, d increases, the Marcenko-Pastur distribution law becomes a 
good approximation of f\(d,j3,z). In particular, the two curves are relatively close for small z, already 
for d = 3. 
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Fig. 1. Comparison between the Marcenko-Pastur distribution and the empirical distribution obtained for j3 = 0.55 and 
d — 1, 2, 3 in the quasi equally space case, and uniform fi(z) 



V. Applications 

Here we present some applications where the results derived in this work can be used. 

The closed form expression of the moments of f\(d,(3,z), given by (1341) , can be a useful basis 
for performing deconvolution operations, as proposed in [8]. As for the asymptotic approximation, we 
show below how to exploit our results for the estimation of the MSE provided by linear reconstruction 
techniques of irregularly sampled signals. 

Let us assume a general linear system model affected by additive noise. For simplicity, consider a 
one-dimensional signal, s(x). When observed over a finite interval, it admits an infinite Fourier series 
expansion [1], [2]. We can think of the largest index M of the non-negligible Fourier coefficients of the 
expansion as the approximate one-sided bandwidth of the signal. We therefore represent s(x) by using 
2M + 1 complex harmonics as 

1 M 

s(x) = ^_ Y a^ lx (19) 
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Now, consider that the signal is observed within one period interval [0,1) and sampled in r points placed 
at positions x = [xo, ■ ■ ■ , x r -i] T , x q £ [0,1), q = 0, . . . ,r — 1. The complex numbers a,£ represent 
amplitudes and phases of the harmonics in s(x). The signal samples s = [s(x ), . . . , s(x r _i)] T can 
be written as s = G^a, where the matrix G is given in (Q]). The signal discrete spectrum is given by 
the 2M + 1 complex vector a = [a-M, ■■■ ,Oq,... ,clm] T ■ We can now write the linear model for a 
measurement sample vector p = [p(xo), . . . ,p(x r -i)] T taken at the sampling points x q 

p = s + n = G f a + n (20) 

where n is a random vector representing measurement noise. The general problem is to reconstruct s or 
a given the noisy measurements p [4], [5]. A commonly used parameter to measure the quality of the 
estimate of the reconstructed signal is the mean square error (MSE). In [l]-[3] it has been shown that, 
when linear reconstruction techniques are used and the sample coordinates are known, the asymptotic 
MSE (i.e., as the number of harmonics and the number of samples tend to infinity while their ratio is 
kept constant) is a function of the asymptotic eigenvalue distribution of the matrix T = /3GG', i.e., 

MSE = E MR (21) 

a |_ASNR m + j3_ 

where the random variable A has distribution f\(d,(3,z) and SNR m is the signal-to-noise ratio on the 
measure. We therefore exploit our asymptotic approximation to f\(d,(3,z) to compute ([2Tb . 

Figure [2] shows the MSE obtained as a function of the signal-to-noise ratio SNR m . The curves with 
markers labeled by "d = 1,2,3" refer to the cases where the signal has dimension d and the sampling 
points are quasi-equally spaced with jitter x, uniformly distributed over [0, 1), and (3 = 0.729. The curve 
labeled by "MP" (thick line) reports the results derived through our asymptotic (d — > oo) approximation 
to the eigenvalue distribution, while the curve labeled by "Equally spaced" (dashed line) represents the 
MSE achieved under a perfect equally spaced sample placement, i.e., when the eigenvalue distribution 
is given by f\(d,/3,z) = 5d{z — 1). Notice that the MSE grows as d increases and tends to the MSE 
obtained by a Marcenko-Pastur eigenvalue distribution. Instead, as expected, the "Equally spaced" curve 
represents a lower bound to the system performance. 

Figure [3] presents similar results but obtained for d = 2 and different values of f3. We observe that the 
MSE obtained through our asymptotic approximation (the curve labeled by "MP") gives excellent results 
for values of ft as small as 0.2, even when compared against the numerical results derived by fixing 
d = 2. For j3 = 0.6 (i.e., when the ratio of the number of signal harmonics to the number of samples 
increases), the approximation becomes slightly looser, and the MSE computed by using the Marcenko- 
Pastur distribution gives an upper limit to the quality of the reconstructed signal. Note that the smaller 
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Fig. 2. MSE as a function of the signal-to-noise ratio for d = 1,2,3. The curves are compared with the results obtained 
through our asymptotic analysis (MP) and with the equally spaced case 



the (3, the higher the oversampling rate relative to the equally spaced minimal sampling rate (3=1. We 
thus observe how our bound becomes tighter as the oversampling rate increases. 

To conclude, we describe some areas in signal processing where the above system model and results 
find application. 

i) Spectral estimation with noise. Spectral estimation from high precision sampling and quantization 
of bandlimited signals uses measurement systems which are usually affected by jitter [20]. In such 
applications the quantization noise corresponds to the measurement noise and the jitter is caused by 
the limited accuracy of the timing circuits. In this case the sampling points are mismatched with 
respect to the nominal values, thus for d = 1 we have: x q = - + ^f- with some sampling rate 1/r. 
Note that the exact positions of the samples are not known and the case studied in this paper (i.e., 
MSE with exact positions) gives a lower bound to the reconstruction error. 

ii) Signal reconstruction in sensor networks. Sensor networks, whose nodes sample a physical field, like 
air temperature, light intensity, pollution levels or rain falls, typically represent an example of quasi- 
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Fig. 3. MSE as a function of the signal-to-noise ratio for (3 = 0.2, 0.6. The curves are obtained for d = 2 and compared 
against both the equally spaced case and the results derived through our asymptotic analysis (MP) 



equally spaced sampling [3], [9], [21], [22]. Indeed, often sensors are not regularly deployed in the 
area of interest due to terrain conditions and deployment practicality and, thus, the physical field is 
not regularly sampled in the space domain. Sensors report the data to a common processing unit (or 
sink node), which is in charge of reconstructing the sensed field, based on the received samples and 
on the knowledge of their coordinates. If the field can be approximated as bandlimited in the space 
domain, then an estimate of the discrete spectrum can be obtained by using linear reconstruction 
techniques [3], [23], even in presence of additive noise. In this case, our approximation allows to 
compute the MSE on the reconstructed field. 
Hi) Stochastic sampling in computer graphics and image processing. Jittered sampling was first examined 
by Balakrishnan in [24], who analyzed it as an undesirable effect in sampling continuous time 
functions. More than twenty years later, Cook [25] realized that the effect of stochastic sampling 
can be advantageous in computer graphics to reduce aliasing artifacts, and considered jittering a 
regular grid as an effective sampling technique. Another example of sampling with jitter was recently 
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proposed in [26], for robust authentication of images. 

VI. Conclusions 

We studied the behavior of the eigenvalue distribution of a class of random matrices, which find large 
application in signal and image processing. In particular, by using asymptotic analysis, we derived a 
closed-form expression for the moments of the eigenvalue distribution. Using these moments, we showed 
that, as the signal dimension goes to infinity, the asymptotic eigenvalue distribution tends to the Marcenko- 
Pastur law. This result allowed us to obtain a simple and accurate bound to the signal reconstruction 
error, which can find application in several fields, such as jittered sampling, sensor networks, computer 
graphics and image processing. 

Appendix A 
Proof of Theorem 13. II 

Using ©, the term TrE [T^] in dD can be written as: 



Tr i i T d) 



X 



E 
x 



E 
x 



E' 

t-x 



-d)u(£ p ),u(lx) 



exp -j2vr^xj (qi) (£i -t 



i+l], 



i=l 



E E 

Le£d QeQd 



exp -j27r^xj (q . ) (£ i -£ [m] ; 



i=l 



where Q d and Cd are sets of integer matrices such that 



{Q | Q = [qi, . . . , q p ], q; = [q i>u 



0,...,p-l} 



C d = {L\L = [£ 1 ,...,£ P \, £ l = [£ i<1 ,...,£ i4 ] T ,e i<m = -M,...M} 



(22) 



and 
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A. Set partitioning 

We now apply the definitions in Section lTlI-Al in order to rewrite (l22l using set partitioning. In particular 
by considering the vector [i = /x(Q) = [/xi, . . . , // P ] T where [ii = /i(qi) and is the i-th column of Q, 
we observe that: 

• the vector /i is uniquely defined by Q, and a given uniquely defines a matrix Q G Qd since //(•) 
is an invertible function; 

• a given /x induces a partition u>(fj,); 

• since r is the number of values that the entries ^ can take, there exist r!/(r — k(fi))\ matrices 
Q G Qd generating a given partition of V made of k{fx) subsets. In other words r!/(r — k(fi))\ 
distinct /n's yield the same partition 

Since the random vectors x^( q ,) and x^( q /,) are independent for q' ^ q", for any given Q the average 
operator in (1221 factorizes into k(fi) terms, i.e., 



E 



exp -j27r^xj (qi) (£i-£| 



i+l] 



i=l 



E 

A' 



exp -j2vr ^ xj. (4 - 



i=l 



fe( M ) 

n« 

fc( M ) 

i=i tj 



exp I -j27rx^. ^ £j - £ [i+1] 



(23) 



indeed, for every i G Vj{n), we have /Uj = 7^. In the last line of d23l . we exploited the following two 
definitions 

C = exp(-j27r/» 

and 



Wo 



» = ^ ii-t]i+i] (24) 
ieP,(/x) 

Also, note that, in the product in (|23l , each factor depends on a single random vector, x 7j . Since 
x /z(q) = q/P + x /i(q)/P an< i m(') i s invertible then, by defining x 7j = /u _1 (7j) we have 



X 7, = ^/P + ^fJP 



and 



E 



(25) 
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In the last term of (1251 ) we removed the subscript jj from the argument of the average operator, since 
the distribution of x 7j . does not depend on jj. Summarizing, the term TrE [T^] in ([8]) can be written as 



TrE[Tg 



1 



e eil 



x ^w 3 -(/i) 



E 

X L 



x T w 3 (fj,) 



(26) 



Since each Q is uniquely identified by a vector /x, we can observe that 

E /(#«) = E E /(/*) = E E E /(/*) (27) 

QeSd oJSfip ^e.M(u>) fc=l we^p.fc /xeA4(a;) 

for every function /(/x). Recall that, in d2"7T ), A^(cj) represents the set of inducing a given partition u>. 

From the definitions in Section UlI-AI it follows that, if /x induces u>, then k(fz) = k(u), Vj[p) = 
Vj(oj), and Wj(fJi) = Wj(u>), j = 1, . . . , k(u). Therefore, 



TVe[T$] 



sEE E EIK 



r p 



r p 



k=l wer2 P , fc (>6M(w) Le£ d j=i 



EE E Enc 

fc=l wef2 Pifc /i£.M(u>) Le£ d j=l 



E 



X T Wj (/x) 



x T w 3 (u>) 



fc=i u>EQ P:k Le£d peA^(oj) 



(a) 1_ 



nc 



x ^,w 3 (w) 



ni[c 



x T w 3 (a;) 



(28) 



fc=i ueti p , t Le£d 



In (|28]) we defined 



fieM(w)-?' =1 



k d 



(29) 



^,L)=n|[c^ M ]=nrif [c 

3=1 j=l m=l 

where x m and Wj m are the m-th entries of x and Wj, respectively. In the equality "(a)" we exploited the 
fact that the term (^"^M does not depend on /x and can be factored from the sum over /jl. As for the 
term X)/zeA4(u>) llj=i C*^™ 3 > we have the following lemma. 

Lemma A.l: Let u G fl Pi k, let wi, . . . , be vectors of size d with integer entries, defined as in (l24b . 
Let M.{u) be the set of vectors /x inducing u>. Then 



e ik w 

/xe7MHi=i 



(30) 



where u(u)') = (— l) fc ' l IIi'=i(l^ ? i'( £ * ; OI ~~ ■"■)'» 7j = "1j(^)> an( ^ wnere ^fc,h is the set of vectors u' of 
size k, representing the partitions of the set V = {1, . . . , k} in h subsets, namely, V[(u'), . . . ,V' h (u'). 
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Proof: The proof can be found in Appendix |B] 
By applying the result of Lemma IA.1I to d28l ), we get 



p k 



T'jM = EE E E 



fc=i fe=i ueJJp,j, cj'efifc.h 



r h u(u') 



j'=i \i'e7V(a>') 



Considering that 



h d 



X[8\ y, **'M =IHlM E 



j'=i \i'e7V(w') 
and by using d29l ) and d3Tl ), we have 



i'=lm=l \i'€P j ,(w') 



(31) 



E^IlM E 



k d h d 



e e nn f 

£ie£i £„eCii=im=i m 

sniMrt'l e *m 

£edj=l j'=l \i'67V(w') 

where the subscript m highlights the dependency of £ on M. 
In conclusion, 

V k h I i\ 

r u(u> ) 



n n 6 [ e ^i'mC") 



Tr | K] = EE E E 



r p 



(32) 



(33) 



fc=i h=i ue!J Pii u/efifc.h 

To compute E[A^ J, we consider the limit in ([8]). By using the definition ((6]), we first notice that 

rP(2M + l) d ~ (2M + i)d(p-M-i) 
Then, by using (|33l in ([8]), we obtain 



p fc 



lim y^y^^ —r, — r—rr u{u}')ibM (.(*>,(*}') 

M,r^+oo ^ ^ (2M + l)d(p-h+l) / YM V 

a fc=i ft=i ue!) Pii o)'en ti/ , 



lim 



E E <«') 

k=l h=l u>eQ P:k ui'&Q k-h 

EEr ft E E «(«x«> 

fe=i /i=i ue!J Pl t u»'ef2 fc ,h 



M^oo (2M + 1)p- /i + 1 



(34) 



The second equality in d34l holds since, for any given p, the sums X^g^ ^ and Yloj'&Q k h 316 over a 



finite number of terms, and the coefficients u(u/) are finite and do not depend on M. Therefore, the 
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limit operator can be swapped with the summations. The coefficient v(u,u') is defined as 

V((jO,U = lim r r—r 

V ' ' M^oo (2M + l)P-h+l 



lim - 

M^oo (2M + 



('-') 



lim 

M^oo (2M + 



; iedj=i j'=i \i'ePy(w) / 

; ted j=i j'=i \i'ep i >(u>) ) 



where, in the equality (a), we introduced the characteristic function of x, defined as C(s) = E[e S2 ]. We 

X 

now consider three possible cases: 

• if h = 1, then fi^j = {[1, . . . , lj}, thus we only consider u' = [1, ... , 1]. Then, Vi(u}') = {1, ... ,k} 

and 



E 

i'eVi((o') 



E 

i'6{l,...,fe} 
k 



i'=l 
k 



E E 



(36) 



i=l 



and by consequence 5 ^X)i'ep./(w') l "i'( a, )l = 1- Hence, 

t,( W ,w') = / []C (-j2vr/? 1 /S'(^)) dy 



(37) 



where, in analogy with (124b . we defined 



w i = E y * ~ y [i+n 

ieVjiw) 

yi e R, i = 1, . . . ,p. We denote by y the vector y = [yx : . . . , y p ] T ; 

if 1 < h < k, the argument of the <5(-) function in 051 ) is always a function of the indices l{. Thus 

\ 

dy 



where 5d(0 denotes the Dirac's delta; 
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if h = k, the cardinality of Q^h = &k,k is S(k, k) = 1 and Qk,k = {[!>■••; &]}■ Thus, we only 
consider uj 1 = [l,...,k]. It follows that: 

„ k k / \ 

v(u,u') = / nc(-j2vr/? 1 /S(u;)) E ffl i'M Uy 

,/Wp i=i i'=i \i'eP 3 -([i,...,fc]) / 

= jf J] C (-j2vr/? 1 /S(u;)) 5 D ( J dy 



k 

' Hp j=l \>'-V,t[l /,.) 

Since Vj([l, ... ,k]) = {j} and C(0) = 1, we have 



(38) 



,(u,,[l,...,k]) = [ f\c(-^(3 1 l d Wj {u))5 D {w 3 {u)) dy 

l[C(0)5 D (wj(u)) dy 
w P i=1 

/ l[S D ( Wj (u)) dy (39) 



As a last remark, if /c = 1, we have h = 1 and = f2 p i = {[!,..., 11}. Then Wj{u) = Y7i=i w i = 0- 



p 



Using (I37T ). we obtain 

[] C (-j2vr/3 1 /^ j (^)) dy = / [] C(0) dy = 1 

. i P j—i J Hp j—i 

Appendix B 
Proof of Lemma [A~T| 

Recall that M.(uj) denotes the set of vectors /i = [/xi, . . . , inducing the same partition us. As defined 
in Section [Tll-Al if u G f2 P) fc, then each /i G M.{u) contains k distinct values, namely, 7 = [71, . . . ,7fc] 
where < jj < r, j = 1, . . . , k and / jj> for each j, j' = l,...,k and j 7^ j'. Therefore, from (IA.1I) 
we can write 

E IK' W E IK w 

/xeX(w)i=i 71 .-".7*^=1 

where the symbol £/yi,...,7* indicates a sum over the variables 71, . . . , 7fc with the constraint that 7j 7^ 7j< 
for every j, / = 1, . . . , k and j ^ j'. Notice that the values jj (j '• = 1, . . . , k) are the scalar counterparts 
of the integer vectors vi, . . . , v^., v.,- = [vji, . . . , Vjd] T , < Vj m < p, m = 1, . . . , d, through the invertible 
function fi(-), i.e., jj = /i(vj), j = 1, . . . , k. Hence, by definition of x, we have x 7j = x M ( Vj ) = vj and 
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in conclusion 

k k 

Y Tl^* j = Y Tic*** = Y e J ™ i+ '" +v ^ k (40) 

We now compute the last term of d40b by summing over one variable at a time. We first notice that, 
for every set vi , . . . , v n of distinct vectors 

r — n w = 



E <T 

V^Vi,...,V„ 



J 3= 

In particular when £ v £ vTw = 0. 

Let us arbitrarily choose the variable v^. If by hypothesis / 0, then by summing (1401 over we 
get 

fe-i 

^VfwH hvjw fc _ _ ^ y ~y s ^vJwiH hv^jWfc-i^vJwfc ($Y) 

Vl,...,Vfc 7=1 Vl,...,Vfc_i 

We compute separately each of the k— 1 contributions in (l4ll . In particular, the generic j'-th term (j = f) 
is given by 

_ hvJ.jWfe-i^v^wj, _ _ ^""^ ^vJwH hvj,(w ;j /+w fc )+vj_ 1 w i ._i 

Vi,...,V fe _i Vi,...,V)b_i 

We now proceed by summing over the variable Vj>. If by hypothesis wy + ^ 0, this summation 
produces k — 2 terms. Again, we consider each term separately. This procedure repeats until a subset S 
of {1, ... , k} is found, such that s = YlieS ™* = ^- 

In this case, the contribution of the n-th sum is given by r — (k — n) where n = \S\ is the cardinality 
of S. Overall, after n sums the total contribution is 

(-iy-i(n-l)!(r-(fc-n)) £ ^ 

v i ,je{i,...,k}-Sje{i,...,k}-S 

The factor (n— 1)! accounts for the number of permutations of the elements in S, once the first element is 
fixed (remember that we arbitrarily chose the first variable of the summation). The factor (— l) n_1 takes 
into account that we summed n—1 times with the condition w / 0, which implies n — 1 sign changes. 
Eventually, the term X)v i je{l,...,fc}-5 IX/ell fc}-.sC v ^ Wj is similar to the last term in (1401 where only 
k — n variables v are involved. 

This procedure repeats until we sum over all variables v. This is equivalent to check if for all possible 
partitions of {!,..., A;} in h subsets V\, ... , Vh, h = 1, . . . , k the condition s± = S2 = • • ■ = Sh = 
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holds, with Sj = Yliev- n i = I^j'I' and X)j n i = ^- m this case, the contribution is given by 

h 

H(-l) n *-\n j -l)\p r (n 1 ,...,n h ) 
j=i 

and it is otherwise. Here p r (ni, . . . ,iih) = (r — (k — ni))(r — (k — n\ — n2)) • • • (? — (k — n\ — rt2 

ra/i-l)). 

In conclusion, we can write 

E r^ i+ - +v ^=E E ti(«0p,(cy) n « ( E 

v i>-»,v* ft=io)'e!] fc , k i'=i \i'e7V(u/) 

where u(c«/) = (-l) fc - h nf'=i(l^7'( w ')l ~ x ) ! and Pr(w ; ) is a polynomial in r of degree h. For large r, 
p r {uj') ~ r ft , thus proving the lemma. 
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